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QUANTUM MOLECULAR DYNAMICS MODELS (QMD) 

The emphasis of this contribution is to high energies, therefore, QMD non-relativistic 
models (E < 2AGeV) are not treated here. 

The first relativistic QMD model was RQMD [1]. This model is no longer supported 
since the year 2000 ,RQMD is used in FLUKA for A-A collisions below 5AGeV. 

A second relativistic model is UrQMD [2], it is used in CORSIKA Cosmic Ray 
cascade code below SOAGeV. 

Let me mention some efforts within the FLUKA collaboration which will not be 
treated here: F.Cerutti et al. add (approximate) energy conservation , evaporation, and 
residual nuclei to RQMD M.V.Garzelli et al. construct a low enery QMD for A-A 
collisions in FLUKA [|4j]. There are efforts in Milano and Houston to construct a fully 
relativistic model for A-A collisions to be inserted into FLUKA. 

A relativistic QMD model is a Lorentz invariant cascade (molecular dynamics) with 
nucleons of both nuclei and all produced hadrons as participants. Properties of such 
models are: (i)A formation zone cascade of all produced hadrons. (ii) Elementary 
interactions used in the models include: (1) h + h — > resonance; resonance + resonance 
and resonance decay (similar to HADRIN in FLUKA), (2) high enery: h + h — > 
hadronic chain; 2 hadronic chains and Lund like chain fragmentation, (3) chain fusion 
(called formation of color ropes) in RQMD, (4) empirical parametrization of all cross 
sections, (5) pQCD description of hard collisions (UrQMD). 



The kinematics of the UrQMD model 

The model is fully described in [2]. It is based on the covariant propagation of all 
hadrons on classical trajectories in combination with stochastic binary scatterings, color 
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string formation and resonance decay. It includes the Monte Carlo solution of a set of 
coupled partial integro-differential equations. Each nucleon is represented by a coherent 
state (/j,c = 1) 



(l>i{x;qi,pi,t)= (^-^^ exp I -^(x- ^,(0)^ + 1^(0^1 (1) 

which is characterized by 6 time-dependent parameters, qi and pi. L, (the extension 
of the wave packet in coordinate space) is fixed. The total n-body wave function is a 
product of coherent states ([1]) ^» = Yli^ii^i^ii Pij)- The Hamiltonian H of the system 
contains a kinetic term and mutual interactions Vjj (H = Y.iTi+ j'Lij^ij)- This yields 
an Euler-Lagrange equation for each parameter. 

d{H) . d{H) 

Pi = -^ and qi = ^. (2) 
oqi dpi 

These are the time evolution equations which are solved numerically. The UrQMD 
Hamiltonian contains : Ekin , , E^lf , E]f '''''' , Ef^"'"'"'' , . 

Please note: As one sees from the Hamiltonian UrQMD is not really a Lorentz in- 
variant molecular dynamics. Therefore, the results of the model might depend strongly 
on the reference frame in which the calculation is done. To minimize this frame depen- 
dence the authors [2] use a frame-independent definition of the cross sections (via using 
the impact parameters in the two-particle rest frames). They give as example for this 
minimization in S-S collisions: The multiplicities and collision numbers vary only by 
less than 3% between the lab and CM frames. 

Let us mention, that RQMD has a manifestly Lorentz invariant eq. of motion. Using 
4-vectors for positions and momenta, each particle carries its own time. The 2N ad- 
ditional degrees of freedom are fixed by 2N constraints: The N mass shell constraints: 
Hi = pj — mf — Vi = 0, (V;-: quasi potential) and (N-1) time fixations. The 2Nth con- 
straint: A relation of times of particles to the evolution parameter T. 

Projectile or target nucleus are modeled according to a Fermi-gas ansatz. The cen- 
troids of the Gaussians are randomly distributed within a sphere with the radius R{A), 
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po is the nuclear matter ground state density. If the phase-space density at the location 

of each nucleon is too high (i.e. the area of the nucleus is already occupied), then 

the location of that nucleon is rejected and a new location is randomly chosen. The 

initial momenta of the nucleons are randomly chosen between and the local Thomas- 

i 

Fermi-momentum: pp"^ = he (Sn^p) \ with p being the corresponding local proton- 
or neutron-density. One disadvantage of this type of initialization: the initialized nuclei 
are not in their ground-state with respect to the Hamiltonian used for the propagation. 
The parameters of the Hamiltonian are tuned to the equation of state of infinite nuclear 
matter and to properties of finite nuclei (such as their binding energy and their root 
mean square radius). One can use a so-called Pauli-potential in the Hamiltonian: This 



has the advantage that the initialized nuclei remain stable whereas in UrQMD with 
the conventional initialization and propagation without the Pauli-potential the nuclei 
start evaporating single nucleons after approximately 20 - 30 fm/c. One drawback 
of this potential: the kinetic momenta of the nucleons are not anymore equivalent to 
their canonic momenta, i.e. the nucleons carry the correct Fermi-momentum, but their 
velocity is zero. The impact parameter of a collision is sampled according to the 
quadratic measure (dW ~ bdb). At a given impact parameter the centers of projectile 
and target are placed along the collision axis in such a manner that a distance between 
surfaces of the projectile and the target is equal to 3 fm. The momenta of the nucleons 
are transformed into the system where the projectile and target have equal velocities 
directed in opposite directions of the axis. After that the time propagation starts. During 
the calculation each particle is checked at the beginning of each time step whether it will 
collide within that time step. 

The relativistic QMD models are compared to a large sample of data in the publica- 
tions of the authors. Up to energies of about ^/s = 200 GeV usually a good agreement 
is found. The experimental collaborations at RHIC compare their data quite often to 
RQMD results. One example: The PHOBOS Collaboration in their white paper [5] find 
a good agreement to dNch / dy from RQMD for Au-Au collisions to their data, the reason 
for this is the chain fusion build from the beginning into RQMD. 

The main problems of the relativistic QMD models, which make them difficult to 
apply for shower simulations are (i) problems with the energy conservation, (ii) the 
missing evaporation and residual nuclei and (iii) the partly excessive computer running 
times. 

DPM AND QGSM MODELS 

The Dual Parton Model DPM and the Quark Gluon String Model QGSM are two 
models, which are largely equivalent in their construction, only with some characteristic 
differences. The DPM is due to Capella, Tran Than Van and collaborators [6.]. The 
most detained Monte Carlo versions of the DPM are PHOJET |0, [sj] for h-h and y-h 
collisions and DPMJET for h-A, A-A and y-A collisions [|90. 

The QGSM is due to Kaidalov, Ter-Martyrosian and collaborators fl^. The 
Monte Carlo version QGSJET for h-h, h-A and A-A collisions is due to Kalmykov, 
Ostapchenko and collaborators [,11,1 . 



The construction of the PHOJET multichain model 

We restrict us in this contribution to describe the PHOJET model, which is used 
directly for h-h collisions and for all elementary Glauber collisions in h-A and A-A 
collisions in DPMJET . There are no essential differences in the formulations of the 
Glauber model between different Monte Carlo models. 

The (soft) Born cross section of the supercritical pomeron has the form Os = 
^2^a(o)-i jj^g supercritical pomeron has a(0) > 1., therefore it clearly violates uni- 



tarity. According to the Froissart bound the cross section asymptotically should not rise 
faster than (log*) . 

If we start to construct the full model, which is unitarized, we should introduce some 
more input Bom cross sections. Very important is the hard cross section, which we 
calculate according to the QCD improved parton model: 



^hard^^^^cutoff) ^ / ^^^^^^^^^ £ 1_ 



f n^\f n^\Zl}±^^£L2lL(A(r. „cutoff^ (a\ 

JaA^l^Q )jbj{X2,Q ) ^ fc)(/7_L-p_L )i (4) 

where fa,i{xi-,Q^) is the distribution of the parton i in a. 

One of the most important difference between PHOJET/DPMJET and QGS JET is in 
the p±_ cutoff p^^toff : PHOJET/DPMJET use p^^^off ^^^^^ ^^^^i energy. QGSJET uses 
^cutoff constant, independent of the energy. 

We introduce furthermore the cross sections for high-mass single and double diffrac- 
tion Od and for high-mass central diffraction Oc according to the standard expressions. 
The amplitudes corresponding to the one-pomeron exchange are unitarized applying 
an eikonal formalism. In impact parameter representation, the eikonalized scattering 
amplitude has the structure 

a{s,B) = i(^Q (l-.-^M)) (5) 

with the eikonal function 

X{s,B) = Xs{s,B)+Xh{s,B) +Xd{s,B) +Xc{s,B). (6) 

Here, Xi{s,B) denotes the contributions from the different Bom graphs: (S) soft part of 
the pomeron and reggeon, (H) hard part of the pomeron (D) triple- and loop-pomeron, 
(C) double-pomeron graphs. 
The eikonals Xii^i^) are defined as follows 

= |Mexp[-^l. (7) 

The free parameters are fixed by a global fit to proton-proton cross sections and elastic 
slope parameters. Once the free parameters are determined, the probabilities for the 
different final state configurations are calculated from the discontinuity of the elastic 
scattering amplitude (optical theorem). 

The total discontinuity can be expressed as a sum of graphs with kc soft pomeron cuts, 
Ic hard pomeron cuts, trie triple- or loop-pomeron cuts, and ric double-pomeron cuts by 
applying the Abramovski-Gribov-Kancheli cutting mles . In impact parameter space one 
gets for the inelastic cross section 

= P|^P|^^» P2£)!lexp[-2,(...)1 (8, 
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with 

/oo 
d^B a{kc,lc,mc,nc,s,B) ^ Gtot (9) 

kc+lc+mc+nc = l 

where Ctot denotes the total cross section 

In the Monte Carlo realization of the model, the different final state configurations 
are sampled from Eq. ([8]). For pomeron cuts involving a hard scattering, the complete 
parton kinematics and flavors/colours are sampled according to the Parton Model. For 
pomeron cuts without hard large momentum transfer, the partonic interpretation of the 
Dual Parton Model is used: mesons are split into a quark-antiquark pair whereas baryons 
are approximated by a quark-diquark pair. The longitudinal momentum fractions of the 
partons are given by Regge asymptotics . We give it here for an event with ns soft and 
"A ("A > 1) hard cut pomerons, sea-quarks are used at the chain ends if we have more 
than one soft pomeron. 

-j^ 2% +2 

p{xi,...,X2n,,-,X2n,+2+nh) ^ ^{ Yl 7)^2^ 

2ns+2+nf, 2ns+2+ni, 

n si^^^QdSii- £ X,). (10) 

!=2n5+3 /=1 

The distributions g{xi, Qi) are the distribution functions of the partons engaged in the 
hard scattering. The momentum fractions of the constituents at the ends of the different 
chains are sampled from this exclusive parton distribution. 

After all this we have all chains defined and PHOJET/DPMJET continues with 
hadronizing all multiple chains using the Lund code JETSET (PYTHIA). 

Now we are able to compare the multichain model Phojet with particle production 
data. There are many comparisons to data published in the PHOJET, DPMJET and 
QGSJET literature. Here we present only two examples: the average multiplicity of all 
kinds of secondary particles in p-p collisions as function of the energy in Figdland the 
rapidity distribution of charged hadrons in central S-S and S-Ag collisions in Figl2l at 
SPS energies. 

Next we present comparisons of the DPM-models (DPMJET-III) with RHIC data 
We first present some comparisons, where DPMJET-III is used in its pre-RHIC form. In 
Figl3]we compare rapidity distributions of charged hadrons in p-p and d-Au collisions 
according to DPMJET with RHIC data from the PHOBOS collaboration. In FigHwe 
compare the transverse momentum distribution measured at RHIC by the PHENIX col- 
laboration with PHOJET calculations, we find the hard collisions very well represented 
in PHOJET. 

For other comparisons DPMJET needs some modifications to get agreement with the 
RHIC data. One of the most important modification is the Percolation of hadronic 
strings in Dpmjet-III 

Using the original Dpmjet-III with enhanced baryon stopping and a centrality of to 
5 % the DPMJET multiplicities are larger than the ones measured in Au-Au collisions 
at RHIC. A new mechanism needed to reduce Nch and dNch/dr}\ri=o in situations with 
a produced very dense hadronic system. We consider only the percolation and fusion of 




FIGURE 1. Average particle multiplicity proton-proton interactions. Phojet results (curves) are com- 
pared to experimental data (symbols). 




FIGURE 2. Rapidity distributions of negative hadrons in central nuclear collisions at 200 GeV/nucleon. 
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FIGURE 3. Pseudorapidity distribution of charged hadrons produced in minimum bias = 200 GeV 
d-Au and p-p collisions. The results of DPMJET are compared to experimental data from the BRAHMS- 
Collaboration and the PHOBOS-Collaboration. At some pseudorapidity values the systematic PHOBOS- 
errors are given. 
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FIGURE 4. Transverse momentum distribution as measured in p-p collisions at ^/s = 200 GeV by the 
PHENIX collaboration at RHIC compared to the calculation by PHOJET 



soft chains (the transverse momenta of both chain ends are below a cut-off = 
2 GeV/c). The condition of percolation is, that the chains overlap in transverse space. 
We calculate the transverse distance of the chains L and K Rl k and allow fusion of 
the chains for Rl k < Rf"^'"" = 0.75 fm. The chains in Dpmjet are fragmented using 
the Lund code. Only the fragmentation of color triplet-antitriplet chains is available in 
Jetset, however fusing two arbitrary chains could result in chains with other colors. 
Therefore, we select only chains for fusion, which again result in triplet-antitriplet 
chains. Examples are: 

(i) A qi — q2 plus a q^—q^ chain become a q\q-i — qiqA chain. 

(ii) A qi — q2q3 plus a <?4 — <?2 chain become a qiq4 — q^ chain. 

(iii) A ^3 — qiq2 plus a q^ — qi plus a qs—qs chain become a ^4 — ^2^5 chain. 

(iv) A q4 — qi plus a q^ — q?, plus a q^ — qi chain become a q^ — q-i chain. 

The expected results of these transformations are a decrease of the number of chains. 
Even when the fused chains have a higher energy than the original chains, the result will 
be a decrease of the hadron multiplicity Nhadmns- In reaction (i) we observe new diquark 
and anti-diquark chain ends. In the fragmentation of these chains we expect baryon- 
antibaryon production anywhere in the rapidity region of the collision. Therefore, (i) 
helps to shift the antibaryon to baryon ratio of the model into the direction as observed 
in the RHIC experiments. 

In Figl5]we compare the pseudorapidity distributions in Au-Au collisions at 200*A 
GeV as measured by the PHOBOS collaboration at RHIC for different centralities with 
the DPMJET-III results obtained with the model including chain percolation and fusion. 

Further RHIC related improvements (not treated here because of the limited space) in 
DPMJET include: 

(i) Modified distributions in PYTHIA: For the fragmentation of soft chains the 
Gaussian transverse momentum distributions in PYTHIA have to be replaced by expo- 
nential ones. 

(ii) Collision scaling in h-A collisions: To obtain collision scaling in DPMJET we 
have to change the sampling of hard chains. 
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FIGURE 5. Pseudorapidity distributions of charged hadrons in Au-Au collisions at y^s)= 200 GeV for 
centralities 0-5 % up to 40-50 %. The points with rather small error bars are from the DPMJET-III Monte 
Carlo with chain fusion as described in the text. The data points are from the PHOBOS Collaboration. 

(iii) Anomalous baryon stopping: New diagrams lead to more baryons in central 
region. 

(iv) Modified diquark fragmentation: Find missing diagram in diquark fragmentation, 
to get Antihyperon to Hyperon ratios into agreement with experiment. 

MODEL COMPARISONS 

This section is based on a talk of D.Heck, Karlsruhe: Comparison of models in 
the CORSIKA Cosmic Ray cascade code at the VIHKOS CORSIKA School 2005, 
Lauterbad, Germany, May 31 - June 5, 2005 and Dieter Heck, private communication. 

High energy models in CORSIKA used for Eiat > 80 GeV include: DPMJET 2.55 
(J.Ranft,Phys.Rev.D51 (1995) 64), NEXUS 2/3 (J.Drescher et al., Phys.Rep.350 (2001) 
93), QGSJET Ol/imil ( S.Ostapchenko, Nucl.Phys.B(Proc.Suppl.)2005), SIBYLL 2.1 
(R.Engel et al.,Proc.26''' ICRC 1(1999)415) . 

Low energy models in CORSIKA used for Ei^t < 80 GeV include: FLUKA 2003 
(only hadron production model) (A.Fasso et al.,Proc. Monte Carlo 2000 (2001)955), 
GHEISHA 2002 (H.Fesefeld, PITHA-85/02 Aachen (1985)), UrQMD 1.3 (S.A.Bass 
et al.,Prog.Part.Nucl.Phys.41 (1998)225). 

In Figl6] we compare the average charged multiplicity in p-p collisions as function 
of the energy between the high energy models. Only the three QGSJET models differ 
strongly from the rest, at high energies the multiplicity increases much stronger than 
in the other models. This is the result of the energy independent cut-off. In FigjT] 
we compare the Xmax of proton and iron induced vertical showers according to the high 
energy models with X^ax data. Such a comparison is hoped to determine finally the 
composition of the highest energy particles in the cosmic radiation. Even, when the 
models differ considerably like in Figl6] in their properties, we find here only modest 
differences in the Xmax predictions. 



600 



400 




10 10^ lo' lO"* 10^ 



FIGURE 6. Table 1 Charged particle average multiplicity in p-p collisions. 
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FIGURE 7. Penetration depth X,„ax for gamma, proton and iron induced vertical showers as function 
of the energy. 



In Table 1 we compare the CPU times of the high energy and low energy models. We 
find only NEXUS and UrQMD to need much more running time than the other models. 



CPU-times (sec) 


for DEC-Alplia lOOOXP 






low energy model 


100000 p-Air coll lOGeV 


high energy model 


10000 p-Air coll 1 PeV 


FLUKA 


181 


DPMJET 2.55 


271 


GHEISHA 2002 


108 


NEXUS 2 


3145 


UrQMD 1.3 


12200 


QGSJET II 


693 






SIBYLL2.1 


186 



SUMMARY AND CONCLUSIONS 



Code comparisons 

* Within 10 years of CORSIKA code comparisons: models have much improved. 



* Accelerator physics oriented code comparisons could help in a similar way. 

* include evaporation particles and residual nuclei. 

* compare also hadron calorimeter performance, produced and residual radioactivity. 
Relativistic QMD models 

* Impressive performance for nucleus-nucleus collisions up to RHIC energies. 

* Missing: exact energy conservation, excited residual nuclei and evaporation, residual 

nuclei. ( Patches to include this into FLUKA). Computer running times of these 
models excessively long. 

* Construct improved relativistic model which includes all properties needed for cas- 

cades at accelerators, this could become a genuine alternative to DPM, QGSM 
models. 

DPM, QGSM models 

* Impressive performance for hadron-hadron, hadron-nucleus, nucleus-nucleus, 

photon-hadron and photon-nucleus collisions up to present collider energies. 

* Improvements through CORSIKA code comparisons. 

* Acceptable agreement of all models up to Auger Cosmic Ray energies. Includes also 

predictions for all cross sections. 

* These are the models which include best evaporation and residual nuclei needed for 

accelerator applications. 
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